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ABSTRACT 

A  plane,  heated  jet  of  a  constant  property  fluid  exiting 
into  a  plane,  straight  channel  of  the  same  fluid  is  analyzed. 
The  effect  of  buoyancy  on  the  velocity  and  temperature  fields 
is  investigated.   Two  jet  orientations  are  considered;  a 
transverse  jet,  which  enters  the  stream  perpendicular  to  the 
main  flow,  and  a  longitudinal  jet,  which  enters  the  stream 
parallel  to  the  main  flow. 

Solutions  for  the  temperature  and  velocity  fields  for  low 
Reynolds  number  flow  are  obtained  using  a  finite  difference 
scheme.   Results  are  presented  for  three  different  values  of 
the  buoyant  force,  holding  other  variables  constant.   For  an 
isothermal  jet,  results  for  the  velocity  field  are  presented 
at  a  higher  Reynolds  number. 
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I.   INTRODUCTION 

Every  commercial  power  generation  plant  releases  some 
waste  heat  to  the  environment.   Water  from  a  nearby  stream, 
lake,  ocean  is  normally  used  as  a  reservoir  for  this 
rejected  heat.   As  the  demand  for  power  increases,  more  and 
more  heat  will  be  discharged  into  these  bodies  of  water. 
The  advent  of  nuclear  power  plants,  which  have  a  lower  ther- 
mal efficiency  than  conventional  fossil-fueled  plants  will 
tend  to  increase  the  amount  of  waste  heat  per  unit  power 
output.   With  the  emergence  of  these  trends,  concern  about 
the  effect  of  these  thermal  discharges  has  begun  to  grow. 

The  utilization  of  a  body  of  water  for  cooling  purposes 
is  normally  carried  out  by  diverting  a  certain  mass  flow 
from  the  body  of  water,  passing  it  through  a  heat  exchanger, 
and  returning  it  at  an  increased  temperature.   Hence,  the 
temperature  of  the  body  of  water,  at  least  in  the  local 
sense,  is  increased. 

From  the  biological  viewpoint,  there  is  no  doubt  that 
raising  the  temperature  of  water  can  be  harmful  to  some  forms 
of  aquatic  life.   For  example,  for  some  fish  a  seasonal 
temperature  rise  of  as  little  as  2°C  can  be  fatal.   Aside 
from  such  obvious  damage,  raising  the  temperature  of  water 
decreases  its  capability  to  retain  dissolved  oxygen,  increases 
the  metabolic  rate  of  fish,  and  affects  their  reproductive 
habits  and  ability  to  resist  disease.   In  addition,  an 
increase  in  environmental  temperature  is  beneficial  to  some 


potentially  harmful  forms  of  plant  life  and  bacteria  [l]. 
Although  the  immediate  effect  of  such  changes  may  not  be 
obvious,  the  normal  equilibrium  of  the  ecological  system 
has  been  disturbed,  and  the  possibility  of  long  range 
damage  does  exist. 

In  order  to  ensure  that  temperature  disturbances  are 
within  safe  limits,  the  precise  nature  of  the  temperature 
disturbance  must  be  known.   The  quantitative  value  of  the 
increase  in  temperature  above  the  undisturbed  temperature 
of  the  given  body  of  water  as  a  function  of  position  and 
time  would  be  of  significant  value  in  analyzing  the  effect 
of  a  thermal  discharge.   Several  studies  related  to  this 
problem  have  been  carried  out. 

The  overall  effect  of  a  thermal  discharge  on  a  given 
body  of  water  is  one  method  of  approaching  the  problem.   A 
method  for  predicting  the  equilibrium  temperature  of  a 
cooling  lake  was  described  by  Sefchovich  [2],   A  study 
relating  plant  parameters  to  the  average  upstream  tempera- 
ture in  a  stream  or  river  was  developed  by  Nahavandi  and 
Campisi  [3] . 

An  experimental  investigation  of  the  temperature  profiles 
in  the  vicinity  of  a  cooling  plant  located  on  Yankee  River 
was  presented  by  Merriam  [4],   Approximate  analytical  models 
were  compared  with  various  experimental  temperature  profiles 
downstream  of  power  plants  situated  on  rivers  by  Polk, 
Benedict,  and  Lahey  [5], 

A  possible  model  of  a  thermal  discharge  is  a  heated  jet 
exiting  into  a  given  environment.   For  the  case  of  thermal 


discharge  into  a  lake  or  ocean  this  model  takes  the  form  of 
a  heated  jet  exiting  into  a  quiescent  fluid.   An  analytical 
solution  for  the  two  dimensional  temperature  and  velocity 
profiles  for  the  case  of  a  heated,  vertical  ,  laminar  jet 
exiting  into  an  infinite  environment  was  developed  by  Brand 
and  Lahey  [6].   Additional  work  on  a  heated  jet  exiting 
into  a  quiescent  fluid  is  now  being  conducted  at  Oregon 
State  University  [7]. 

The  case  of  a  heated  discharge  into  a  stream  or  river 
may  be  modeled  by  a  heated  jet  exiting  into  a  moving  stream. 
An  experimental  investigation  of  a  buoyant,  turbulent, 
circular  air  jet  exiting  into  a  cross  wind  was  conducted  by 
Ramsey  and  Goldstein  [8].   An  analytical  study  of  the 
surface  spread  of  a  submerged,  buoyant  incompressible  jet 
was  made  by  Tulin  and  Schwartz  [9], 

The  problem  undertaken  in  this  study  was  that  of  a 
heated  jet  in  a  stream.   This  particular  problem  can  be 
approached  in  several  ways.   The  general  problem  is  a  tran- 
sient problem  in  which  the  temperature  and  velocity  fields 
vary  in  three  dimensions.   Factors  which  must  be  considered 
are  diurnal  temperature  and  heat  flux  variations,  tidal 
effects,  seasonal  temperature  and  flow  variations,  atmos- 
pheric conditions,  stream  meander  path,  bottom  contour 
variations,  fluid  property  variations,  the  possible  turbulent 
nature  of  the  effluent,  and  temperature  stratification 
effects.   Because  of  the  complexity  of  this  problem,  several 
simplifications  were  made. 
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The  transient  conditions  were  assumed  to  be  "slow" 
compared  to  the  response  time  of  the  stream.   Therefore, 
at  any  instant  in  time  the  stream  could  be  considered  to 
be  at  a  quasi-equilibrium  state.   The  time  dependent  terms 
were  therefore  not  considered.   The  fluid  properties  were 
assumed  to  be  constant,  for  simplification,  and  the  axis  of 
the  river  was  assumed  to  be  straight. 

It  was  decided  to  attempt  to  isolate  the  effects  of 
buoyancy  on  the  temperature  and  flow  patterns.   Also,  to 
put  the  problem  in  a  more  tractable  form,  it  was  decided  to 
consider  variations  in  two  dimensions  only.   Assuming  that 
one  of  the  dimensions  would  be  the  downstream  dimension, 
there  are  two  possible  choices:   a  plane  parallel  to  the 
river  surface  or  a  vertical  plane  passing  through  the  axis 
of  the  river.   Since  the  buoyancy  forces  act  in  the  vertical 
plane,  the  vertical  plane  was  chosen.   In  order  to  isolate 
buoyant  effects,  the  bottom  of  the  channel  was  assumed  to  be 
straight  and  parallel  to  the  river  surface.   Also,  for  sim- 
plification, the  river  and  the  jet  were  assumed  to  be 
laminar  in  nature.   The  desired  solutions  were  the  tempera- 
ture and  velocity  profiles  as  a  function  of  position  and 
significant  non-dimensional  parameters. 
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II.       ANALYSIS 

A.  COORDINATE  SYSTEM 

The  problem  under  consideration  in  this  work  is  flow 
in  a  channel  into  which  a  known  quantity  of  fluid  is  injected. 
The  injected  fluid  is  at  a  higher  temperature  than  the  fluid 
in  the  channel.   Since  it  is  desired  to  study  the  effects  of 
buoyancy,  variations  in  the  vertical  plane  only  are  consid- 
ered.  The  axis  of  the  river  is  assumed  to  be  straight. 
Rectangular  cartesian  coordinates  are  used,  as  shown  in 
Figure  1. 

Two  different  orientations  of  the  injected  mass  with 
respect  to  the  main  flow  are  considered.   The  longitudinal 
jet  is  injected  parallel  to  the  direction  of  main  flow. 
The  transverse  jet  is  injected  perpendicular  to  the  direction 
of  main  flow.   For  the  longitudinal  jet,  the  origin  of  the 
coordinate  system  is  located  at  the  channel  bottom,  immedi- 
ately below  the  jet  location.   For  the  transverse  jet,  the 
origin  of  the  coordinate  system  is  located  far  enough 
upstream  of  the  jet  so  that  the  incoming  temperature  and 
velocity  profiles  are  not  disturbed  (see  Figure  1) . 

B.  GOVERNING  EQUATIONS 

Applying  the  Boussinesq  approximation,  which  assumes 
that  all  density  variations  other  than  those  giving  rise  to 
the  buoyant  forces  may  be  neglected,  and  neglecting  viscous 
dissipation,  the  steady  two  dimensional  forms  of  the  conti- 
nuity, Navier-Stokes ,  and  energy  equations  are: 
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is.  +  Ix   =    o  (i) 

3x         3y 

2  2 

r       3u         ..    3ut                    3p    ,        r3    u    ,     3    u-i  /ox 

p[u   —  +   v   — J      = &  +    uL 5-  +   oJ  (2) 

3x  3y  3x  3x^         3y^ 
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pc[u  II  +  v  II]      =      k[l!|  +  if|]  (4) 

3x      3y         3x2    3y^ 


where  TQ  is  the  incoming  stream  temperature. 

By  cross-differentiating  with  respect  to  x  and  y,  and 
by  subtracting  equation  (3)  from  equation  (2),  the  pressure 
p  is  eliminated.   The  resulting  equation  is: 

ui?H_+vi5i-ul^-vl^  =   »[4-  +  ^5  -  d-h  -  2^j-ge  S 

3x3y  3y^  3x2  3x3y  3x23y       3yJ        3x^        3x3y2  3x 

(5) 
The   continuity   equation  may   be    eliminated   by   defining 
a   stream   function   i> ,    such   that 

o    -   14       ,    v   *    -  il 

3y  3X 

Only  two  unknown  quantities  remain  to  be  determined, 
i>   and  T.   The  two  remaining  equations  are 


3^   r3V  24>      di\>  r3V2<K        4        3T              .,, 

_l   r L] 1    [ 1]   =   VV  ili  -  g6  —  (6) 

3y   L  3x  J     3x  L  3y  J         ■  *    y   3x 

li  II  -  ii  21  =  a  V2T                                  (7) 

3y  3x    3x  3y 
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Equation  (6)  is  the  momentum  equation  expressed  as  the 
transport  of  vorticity,  where  the  vorticity  is  defined  as 

The  equations  will  be  non-dimensionalized  by  introducing 
the  following  variables: 


U   = 


-   u 


Uo 


v  =  Y_ 

V. 


X   =   * 


D 


Y      = 

D 

0       = 

T    -    TQ 

Tj    "    T0 

$*    = 

U0D 

where  U   is  the  average  undisturbed  velocity,  D  is  the 
channel  depth,  and  T^  is  the  jet  temperature. 

By  substituting  these  into  equations  (6)  and  (7) , 


dill*  3V1*        3^*      3  7  i|;*        1    4      Gr   3  6       /ON 

_I_  (- L-]   -  _1_   [ T_]   =  —  v%*  -  — -  _     (8) 

3Y      3X         3X       3Y         Re         Rez    3X 


3^*  39  _  3if;*  98   _   1_   y2e  (9) 

3Y   3X    3X   3Y      RePr 
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where  Re  is  the  Reynolds  number  based  on  the  channel  depth, 
Pr  is  the  Prandtl  number  of  the  fluid,  and  Gr  is  the  Grashof 
number  based  on  the  channel  depth  and  the  temperature  differ- 
ence between  the  jet  temperature  and  the  incoming  stream 
temperature . 

C.   BOUNDARY  CONDITIONS 

Examination  of  the  governing  equations  reveals  that 
twelve  boundary  conditions  are  necessary  to  completely 
specify  the  problem.   The  necessary  conditions  are  four 
conditions  for  ty   on  x,  four  for  |i  on  y,  two  for  T  on  x  and 
two  for  T  on  y.   These  conditions  will  be  specified  on  the 
channel  bottom,  the  river  surface,  upstream  and  far  down- 
stream of  the  jet.   The  boundary  conditions  for  \p   are  first 
expressed  in  terms  of  velocity,  and  then  related  to  the 
stream  function  from  the  definition  of  \p ,  i.e.,  QA   =  u, 

— -   =  -v. 
3x 

At  the  upstream  boundary,  the  undisturbed  incoming 
horizontal  velocity  profile  is  assumed  to  be  a  20%  truncated 
parabola,  which  is  a  common  assumption  in  stream  flow  [10]. 
A  20%  truncated  parabolic  flow  profile  has  the  maximum 
velocity  located  two  tenths  of  the  total  depth  below  the 
free  surface  (see  Figure  2) .   For  the  longitudinal  jet, 
this  profile  is  modified  by  assuming  that  the  stream 
velocity  is  equal  to  the  jet  velocity  across  the  jet  width 
(see  Figure  3) .   For  the  transverse  jet,  the  upstream 
boundary  is  located  far  enough  upstream  so  that  the  velocity 
profile  is  equal  to  the  undisturbed  velocity  profile.   The 
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Figure  2  -  20%  truncated  parabolic  velocity  profile. 
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Figure  3  -  Velocity  profile  at  upstream  boundary 
for  longitudinal  jet. 
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vertical  component  of  the  stream  velocity  is  assumed  to  be 
zero.   For  the  longitudinal  jet,  the  stream  temperature  is 
assumed  to  be  equal  to  T0  except  across  the  jet  width, 
where  the  temperature  is  equal  to  the  jet  temperature,  T^ . 
For  the  transverse  jet,  the  upstream  boundary  is  located 
far  enough  upstream  so  that  the  incoming  temperature 
profile  is  also  undisturbed. 

On  the  channel  bottom,  the  no  slip  boundary  condition 
requires  that  u  =  0.   The  vertical  component  of  velocity 
is  also  equal  to  zero,  except  at  the  transverse  jet.   In 
that  case,  v  is  equal  to  zero  along  the  river  bottom  except 
across  the  jet  width,  where  v  is  equal  to  the  jet  velocity. 
The  river  bottom  is  assumed  to  be  adiabatic,  yielding  a 
temperature  gradient  boundary  condition. 

On  the  channel  surface,  it  is  assumed  that  there  are 
no  waves;  i.e.,  v  =  0.   By  assuming  that  the  pressure  is 
constant,  and  that  the  river  surface  is  level,  the  appli- 
cation of  the  Bernoulli  equation  along  the  channel  surface 
leads  to  the  result  that  the  surface  velocity  Us  is 

constant.   Since  U   is  constant,  and  the  value  of  U-  at 

s  a 

the  upstream  boundary  is  known,  the  surface  velocity  is 
equal  to  the  undisturbed  surface  velocity.   At  the  surface, 
the  conductive,  convective,  and  evaporative  nodes  of  heat 
transfer  to  the  atmosphere  are  lumped  into  an  overall  heat 
transfer  coefficient.   It  is  assumed  that  radiative  heat 
fluxes  can  be  modeled  by  a  known  constant  heat  flux  q. 
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A  simple  heat  balance  shows  that  the  heat  conducted  into 
the  surface  plus  the  radiative  heat  flux  into  the  surface 
equals  the  heat  convected  out. 

The  velocity  profile  at  the  downstream  boundary  is 
assumed  to  be  parabolic  in  form  such  that  the  total  mass 
flow  is  equal  to  the  undisturbed  mass  flow  plus  the  mass 
flow  of  the  jet.   The  vertical  component  of  velocity  is 
assumed  to  be  equal  to  zero.   It  is  also  assumed  that  all 
heat  added  by  the  jet  is  convected  out  of  the  channel 
through  the  surface.   Hence,  the  temperature  profile 
returns  to  its  original  constant  value  T  .   The  downstream 
boundary  is  situated  far  enough  downstream  so  that  the 
above  conditions  are  met  within  a  given  tolerance. 

The  mathematical  forms  of  these  conditions  are: 

3^L  Uo_  .„.    ,„.2. 


x=0     uL  =  -±  =  u±(y)  =  -£[8(£)-  5(1)  ]        oiy<y-iL 
h      ay  3    d    d  lh 

YjL<y<yju 


-  u  . 
3 


=  ui(y)  yjU<yiD 


uT  =  =  u^  (y) 

ay 


v  =  -  _1  =  0 
3x 


TL  "  To  0<y<yjL 

-  Tj  yjL<y<yju 


=  T 


o  yju<y<D 
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1    o 
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where  the  L  subscript  refers  to  the  longitudinal  jet,  the 
T  subscript  refers  to  the  transverse  jet,  the  i  subscript 
refers  to  the  incoming  velocity  profile,  the  D  subscript 
refers  to  the  velocity  profile  at  the  downstream  boundary, 
y^L  is  the  y-location  of  the  lower  edge  of  the  jet,  y-ju 
is  the  y-location  of  the  upper  edge  of  the  jet,  x-t  and 
x.y  are  defined  similarly  for  the  x-location  of  the  jet, 
and  T   is  the  ambient  air  temperature. 

Since  the  problem  is  cast  in  terms  of  the  stream 
function  ty ,    it  is  desired  to  reduce  the  boundary  conditions 
to  a  more  convenient  form.   At  the  undisturbed  upstream 
boundary,  the  20%  truncated  parabola  gives  the  following 
form  for  the  velocity  profile. 


u  =  11  =   Ha  [s  z  -  s  (i)2] 

8y     3     D      D 


Integrating  from  0  to  D,  it  was  found  that  the  average 
velocity  UQ  is  given  by 

Now,  integrating  u  (y)  from  zero  to  y,  assuming 
*(x,0)  =  0 

i>        =   _°_  [12  (I)2  -  5  (X)3] 

This  is  the  upstream  boundary  condition  for  #  for  the 
transverse  jet.   For  the  longitudinal  jet,  if  it  is  assumed 
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that  the  jet  diameter  is  small  with  respect  to  the  river 
depth,  the  expression  for  \\>   becomes: 

*  =  *u(y)  °iy<YjL 

=  *u(y)  +  UjD^LIiL)        yjL<y<yju 
=  *u(y)  Yju<yiD 

where  D^  is  the  jet  width,  y.j.  is  the  location  of  the  lower 
edge  of  the  jet,  and  y.*,  is  the  location  of  the  upper  edge 
of  the  jet. 

For  the  longitudinal  jet,  the  boundary  condition  at 
the  channel  bottom  may  be  satisfied  by  assuming  41  is  con- 
stant.  Hence,  \p    is  chosen  to  be  zero  on  the  channel 
bottom  for  this  case. 

For  the  transverse  jet,  the  boundary  condition  for  v 
may  be  integrated,  giving 


^T  =  0  for     x<xjL 


x-x  • 


i>T   =  -U.D.  ( 3L) 


X -t  <x<x. 


.u.     ^—^->  *jL^*jU 


tf>T  =  "UjDj  x>xjU 

The  boundary  condition  for  v  at  the  surface  may  also 
be  satisfied  by  choosing  ip  =  constant.   The  value  of  the 
constant  is  equal  to  the  value  of  \p   at  the  surface  as 
given  by  the  upstream  boundary  condition. 

The  downstream  boundary  condition  for  t|i  may  be  found  by 
integrating  the  downstream  velocity  profile,  and  evaluating 
constants  from  the  surface  and  bottom  values  for  ty . 
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After  performing  the  indicated  operations,  and  non- 
dimensionalizing,  the  final  form  of  the  boundary  conditions 
becomes 


X=0      *T*(Y)  =  *U*(Y)  =  ±    (12Y2  -  5YJ) 


^U*(Y) 


Y-Y 


°iY<YjL 


^L*(Y)  "  <  ^U*(YjL)  +  UjDj*{ —]       YjL<Y<YjU 


D--* 


i|>„*(Y)  +  U^D.,* 


'U 


3] 


Yju<Y<l 


3^* 
9X 


=   0 


T 


=   0 


0 

0<Y<Y.T 

—  —  JLi 

1 

YjLlYlYjU 

0 

Yju<Y<l 

Y=0 


<J>L*   =   ° 


ipT 


X-X. 
-U-D.*( Hi) 

1    3  DJ 

v-ujDj* 


°lx<XjL 


XjL<X<XjU 


Xju<x 


3^* 
3Y~ 


=   0 


36 

Ty" 


=   0 
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Y=l 


86T 

=       0  0<X<X- 


3Y 


-     "  JL 


eT       =     l  xjL<x<xju 

=     0  x.„<x 


3Y 

34>*      =      9_ 
3Y  7 


*T*      =      1 


^L*      =      1    +    UjD. 


* 


I!      _    hD         =       hD     (Q    _  } 

3Y  K  K  A 


X+°°  ^L*       =      ^D*(Y)       =       (12    +    SUjD^JY2    +     (2UjDj*    -    5)  Y3 


*T*      =      ^D*(Y)     ~   UjDj 


IT  "    ° 

e         =o 

where  the  L  and  T  subscripts  have  the  same  meaning  as 
defined  previously. 
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III.   METHOD  OF  SOLUTION 

A.   GOVERNING  EQUATIONS 

The  governing  equations  for  the  problem  under  considera- 
tion are: 

di>*    jvV    3^*  fdv2j>*  1   4     Gr  ae 

( )  -  ( )    =   V^*  -  7?   (10) 

3Y     3X       3X      3Y         Re         Rez    3X 


3**  30     3U;*  36        It 

-JE 2 =   v^e  (11) 

3Y   3X    3X   3Y      RePr 


In  order  to  solve  this  problem  using  a  finite  differ- 
ence method,  a  finite  difference  operator  for  each 
differential  operator  in  the  governing  equation  is  required. 
The  method  used  to  determine  these  operators  is  basically 
the  same  as  that  used  by  Dias  [ll].   The  operators  developed 
in  this  study,  however,  are  applicable  to  rectangular  grid 
networks  as  well  as  to  square  grid  networks. 

The  procedure  for  determining  the  finite  difference 
operators  is  to  first  define  an  interpolating  polynomial 
P  (x,  y)  having  n  undetermined  coefficients.   This  poly- 
nomial can  be  used  as  an  approximate  to  any  general  function 
F.   (F  can  be  interpreted  as  a  stream  function,  for  example.) 
By  requiring  that  P  (x,  y)  be  equal  to  F  at  each  of  n 
points,  n  equations  involving  the  n  undetermined  coeffi- 
cients are  generated.   Hence,  the  coefficients  of  P  (x,  y) 
can  be  determined  in  terms  of  the  value  of  F  at  each  node. 
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By  applying  the  desired  differential  operator  to  P  (x,  y) , 
an  approximation  for  the  operator  in  terms  of  the  value  of 
F  at  each  node  is  obtained. 

To  develop  the  biharmonic  operator,  a  polynomial  must 
be  chosen  such  that  the  fourth  partial  derivatives  with 
respect  to  both  x  and  y  are  non-zero.   The  grid  network 
through  which  P  (x,  y)  is  to  be  passed  is  shown  in  Figure 
4.   Since  there  are  25  points  in  this  grid,  25  unknown 
coefficients  are  required  in  the  polynomial  P  (x,  y) . 
Hence,  an  appropriate  choice  for  P  (x,  y)  is 

P(x,y)  =  C]^y4x4  +  C2y3x4  +  C3y2x4  +  C4yx4  +  C5x4 
+  C6y4x3  +  C7y3x3  +  C8y2x3  +  C9yx3  +C10x3 
+C11y4x2  +Ci2y3x2  +C13y2x2  +C14yx2  +C,rX2 

+c16y4x   +c17Y3x   +Ci8Y2x   +C19yx   +C20x 

+C21y4    +C22Y3    +C23y2    +C24y    +C25         (12) 

This  equation  may  be  written  in  matrix  form  as: 

P(x,y)   =   <Vi(x,y)>    {Ci}        i=l,2,...,25         (13) 

where  V.  (x,  y)  is  x^y^  corresponding  to  C^  in  equation  (12) 

This  polynomial  is  now  required  to  pass  through  the 
value  of  F •  at  each  node  j ,  where  F •  is  the  value  of  the 
function  F  at  node  j .   At  each  node  j ,  the  following 
relation  is  obtained: 

P(Xj,yj)  =  <Vi(Xj,yj)>    {C±}        =   Fj      i  =  1,2,.. .,25 

(14) 
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By  applying  this  relation  at  each  of  the  25  grid  points, 
the  following  set  of  equations  results: 


25 


<vi(x1,yi)> 
<Vi(x2  y2)> 


<Vi(x25,y25)> 

-4 


(15) 


-L  f  £-  f    «  •  •  f   Z.  Z) 


Or,  in  matrix  notation: 


{F.}   =   [M]{Cj} 


J       J-/'-f«»»/^--J 


(16) 


Solving  this  matrix  equation  for  C  •  , 


{Cj}   =   [M]"1{Fj} 


_J   """"   X^^»^»»«/^3 


The  coefficients  C_j  of  P  (x,  y)  are  now  known,  and 
P  (x,  y)  can  be  used  as  an  approximation  for  F.   The 
biharmonic  operator  may  now  be  applied  to  the  polynomial 
P  (x,  y) . 

V4(P(x,y))   =  <V4  (  Vi(x/y)>  {C±}) 

v4<(Vi(x,y))>  {Ci}     i  =  1,2, ...,25 


Since  an  operator  for  the  biharmonic  operator  for  the 
central  node  is  desired,  this  expression  is  evaluated  at 
x13,y13. 

V4(P(x,y))       =      <V4     (Vi(x,y))>    [m]"1    {Fi>       i=l,2,...,25 


=      <K 


i>       (F-) 


X~~  X  ^  ^-  /  •   •  •   /^O 


(17) 
(18) 
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where  <K^>  is  a  row  vector  of  constants  obtained  by 
multiplying  the  matrix  expressions  appearing  in  equation 
(17) .   The  result  is  that  K^  are  the  coefficients  of  a 
finite  difference  approximation  for  the  biharmonic  operator 
applied  at  the  central  node  of  the  grid  in  Figure  4. 

In  order  to  ensure  geometric  consistency/  this  same 
interpolating  polynomial  was  used  to  develop  the  operators 
appearing  on  the  left  hand  side  of  equation  (10) .   The 
desired  operator  is  applied  to  the  polynomial,  and  evalu- 
ated at  the  central  node.   The  results  for  a  =  5  are 
presented  in  Figures  5,  6,  and  7,  where  a  is  the  ratio  of 
the  spacing  in  the  x-direction  to  the  spacing  in  the 
y-direction. 

A  finite  difference  approximation  for  the  harmonic 
operator,  which  is  required  in  equation  (11)  is  available 
in  the  literature  [12],   Consistent  approximation  for  the 
other  operators  appearing  in  this  equation  are  also 
available.   These  operators  are  presented  in  Figures  8 
and  9. 

B.   BOUNDARY  CONDITIONS  IN  FINITE  DIFFERENCE  FORM 

Assuming  there  is  no  jet  present,  the  boundary  condi- 
tions for  ty*   are  all  of  the  form 

**   =  fit) 

"St  =  g(t) 

an 

where  t  is  the  space  coordinate  tangent  to  the  boundary, 
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Figure  8  -  Finite  difference  approximation  for 

V2  operator  applicable  at  central  node 
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n  is  the  coordinate  normal  to  the  boundary,  and  f  (t)  and 
g  (t)  are  the  values  specified  at  each  particular  boundary. 
The  first  condition  on  ty*    is  satisfied  in  the  finite 
difference  scheme  by  evaluating  f  (t)  at  each  node  on  the 
given  boundary.   These  values  are  then  used  in  the  finite 
difference  computations. 

To  satisfy  the  boundary  condition  for  the  normal  deriva- 
tive of  ty ,  g  (t)  is  evaluated  at  each  node  on  the  boundary. 
A  finite  difference  approximation  for  the  normal  derivative 
is  then  required.   An  offset  three  point  operator  is  used 
because  the  truncation  error  is  consistent  with  the  trunca- 
tion error  of  the  governing  equation.   The  finite  difference 
operators  used  at  the  boundaries  are  presented  in  Figures 
10  and  11. 

The  temperature  boundary  conditions  assume  various 
forms.   At  the  upstream  and  downstream  boundaries,  the  value 
of  the  temperature  at  the  boundary  is  specified.   At  nodes 
located  on  these  boundaries,  the  known  value  of  the  temper- 
ature is  used  in  the  finite  difference  computations.   On 
the  channel  bottom,  the  derivative  of  the  temperature  is 
specified.   At  nodes  on  this  boundary,  the  offset  three 
point  operator  utilized  previously  is  applied.   At  the 
channel  surface,  the  sum  of  the  temperature  and  its  deriva- 
tive is  specified.   At  nodes  on  the  channel  surface,  a 
mixed  boundary  condition  exists.   In  this  case,  the  value 
of  the  temperature  at  the  node  located  on  the  surface  is 
added  to  the  three  point  operator  shown  in  Figure  11. 
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For  most  computations,  the  jet  center  is  located  between 
two  nodes.   In  addition,  the  jet  diameter  is  assumed  to  be 
less  than  the  grid  spacing  in  all  cases.   The  effective 
diameter  of  the  jet  is  equal  to  the  grid  spacing  on  the 
boundary  on  which  the  jet  is  located.   This  is  a  result  of 
the  fact  that  any  length  smaller  than  the  grid  spacing 
cannot  be  accurately  modeled  in  a  finite  difference  network. 
This  can  best  be  illustrated  by  example.   Figure  12  depicts 
two  possible  jet  configurations,  both  having  the  same  volume 
flow  rate.   Also  shown  are  the  stream  function  profiles, 
superimposed  on  a  grid  network.   In  both  cases,  the  jet 
width  is  less  than  the  grid  spacing.   It  can  be  readily 
seen  that  the  values  of  the  stream  function  at  the  two 
nodes  on  either  side  of  the  jet  are  the  same,  regardless 
of  the  jet  diameter.   The  value  of  the  stream  functin  is 
changed  by  the  amount  of  the  volume  flow  rate  of  the  jet. 
If  the  jet  width  is  smaller  than  the  grid  spacing,  the 
effect  of  the  jet  on  the  finite  difference  grid  is  a  fixed 
flow  rate  between  two  grid  points.   The  effective  diameter 
of  the  jet  is  then  equal  to  the  grid  spacing,  and  the 
effective  jet  velocity  is  the  flow  rate  divided  by  the  grid 
spacing. 

By  modeling  the  jet  in  this  manner,  the  edges  of  the 
jet  fall  on  nodes  of  the  finite  difference  grid.   In  a 
continuous  system,  there  is  a  step  change  in  both  the 
temperature  profile  and  the  velocity  profile  at  these 
points.   This  step  change  cannot  be  modeled  precisely  in  a 
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Figure  12  -  Velocity  profiles  and  resulting  stream 
function  profiles  for  two  jets  located 
between  nodes  1  and  2. 
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finite  difference  scheme.   By  specifying  the  temperature 
at  these  two  nodes  in  the  finite  difference  grid,  the 
temperature  profile  is  as  shown  in  Figure  13.   It  can  be 
seen  that  the  total  area  under  the  temperature  profile 
curve  is  equal  to  the  area  under  a  step  change  of  tempera- 
ture 2  x  Ts .   The  energy  input  of  the  profile  shown  in 
Figure  13  is  equal  to  that  of  the  step  change  of  height 
2  x  Ts.   This  temperature  may  be  defined  as  the  effective 
jet  temperature. 

If  the  three  point  finite  difference  operator  shown  in 
Figure  9  is  used  to  calculate  the  velocity  profile  based 

on  the  relation  V  =  2^L_  ,  it  can  be  shown  that  the  effec- 

3Y 

tive  jet  velocity  defined  earlier  is  consistent  with  the 
definition  of  the  effective  jet  temperature. 

The  above  examples  are  directly  applicable  to  the 
transverse  jet.   By  analogy,  the  results  may  easily  be 
extended  to  the  longitudinal  jet. 

For  the  results  used  to  construct  the  convergence  plot, 
the  jet  center  is  located  at  a  node.  For  this  orientation, 
the  effective  jet  diameter  is  twice  the  grid  spacing. 

C.   NUMERICAL  SOLUTION  SCHEME 

In  the  numerical  solution  of  this  problem,  the  governing 
equations  are  rearranged 

y4^*   =   _  £r  3_e_  +  j-3jj/*  3(V2^*)  _  ^*_   3  ( vV)  ]Po  (18) 

Re  3X     3Y     3X       .3X     3Y 


V26        =      RePr    [li_  £JL  -    d^      £JL]  (19) 

3Y       3X  3X       3Y 
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Figure  13  -  Temperature  profile  over  finite  difference 
grid  in  the  vicinity  of  the  transverse  jet. 
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These  equations  may  be  expressed  in  an  alternate  form  if 
we  recognize  that  the  Jacobian  operator  is  defined  as 


T/,  _*       3A  3B    3A  3B 
J(A,B)   = 

3X  3Y    3Y  3X 


Equations  (1)  and  (2)  may  then  be  rewritten  as 


V%*   =   -  ||  ||  -  JU*,vV)Re  (20) 


72e   =  RePr  J(\\j*,e)  (21) 


In  the  solution  of  these  equations,  an  iterative  scheme 
based  on  the  uncoupled  linear  homogeneous  forms  of  these 
equations  is  used.   Once  this  linear  system  of  equations  is 
solved,  the  right  hand  side  of  equations  (20)  and  (21)  may  be 
calculated. 

The  linear  equations  are  then  resolved,  using  the  calcu- 
lated values  of  the  right  hand  side  of  the  equations  as  known 
inhomogeneities .   This  procedure  is  repeated  until  the  desired 
convergence  criteria  is  satisfied. 

The  uncoupled  linear  homogeneous  forms  of  equations  (20) 
and  (21)  are: 

V  \*      =   0 

v4e   =  o 
The  method  of  solving  these  equations  with  the  boundary 
conditions  for  the  transverse  will  be  illustrated  using  the 
grid  in  Figure  14  as  an  example. 
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The  boundary  conditions  on  ty*    specify  the  value  of  4>*   on 
all  boundaries.   Therefore,  the  value  of  i>*   at  nodes  located 
on  the  boundaries  is  known.   The  value  of  ty*   at  each  interior 
node  is  unknown.   In  order  to  solve  for  the  value  of  ty*   at 
interior  nodes,  one  equation  for  each  node  is  required. 

For  all  ty*   values  not  adjacent  to  any  boundary  (nodes  15, 
16,  21,  and  22  in  example),  an  equation  is  obtained  by  applying 
the  governing  equation  at  the  node.   The  finite  difference 
operator  shown  in  Figure  5  is  used. 

For  all  i>*   values  located  at  nodes  adjacent  to  the  upstream 
boundary  (nodes  8-11  in  example) ,  an  equation  is  obtained  by 
applying  the  boundary  condition  -r^—   =0   at  the  adjacent 
boundary  node  (nodes  2-5  in  example) .   The  finite  difference 
operator  shown  in  Figure  10  is  used. 

For  \l>*   values  at  nodes  adjacent  to  the  downstream  boundary 
(nodes  26-29  in  example) ,  an  equation  is  obtained  by  applying 
the  boundary  condition  ^—-   =0   at  the  adjacent  boundary 
node  (nodes  32-35  in  example) .   The  finite  difference  operator 
shown  in  Figure  11  is  used. 

For  \p*   values  at  nodes  adjacent  to  the  channel  bottom,  but 
not  adjacent  to  either  the  upstream  or  downstream  boundary 
(nodes  12,  20  in  example),  an  equation  is  obtained  by  applying 
the  boundary  condition  t-^—  =0   at  the  adjacent  boundary  node. 
The  finite  difference  operator  shown  in  Figure  10  is  used. 

For  >|;*  values  at  nodes  adjacent  to  the  channel  surface,  but 
not  adjacent  to  either  the  upstream  or  downstream  boundary, 
(nodes  17,  23  in  example),  an  equation  is  obtained  by  applying 
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3  il>*    9 
the  boundary  condition  -^5—  -  *  at  the  adjacent  boundary 

ox       / 

node  (nodes  13,  19  in  example).   The  finite  difference  operator 
shown  in  Figure  11  is  used. 

Since  the  value  of  ip*  is  known  at  nodes  located  on  the 
boundary,  whenever  a  finite  difference  operator  involves  the 
value  of  ty*   on  a  boundary,  the  known  value  of  4>*    is  used. 
The  resulting  system  of  equations  is  a  set  of  algebraic  equa- 
tions for  each  unknown.   This  system  of  equations  is  solved 
by  Gaussian  elimination. 

The  boundary  conditions  for  8  specify  the  value  of  6  at 
the  upstream  and  downstream  boundaries,  and  at  the  two  nodes 
on  either  side  of  the  jet.   However,  the  value  of  the  tempera- 
ture on  the  channel  surface  is  unknown,  and  the  value  of  the 
temperature  on  the  channel  bottom,  except  for  the  nodes  on 
either  side  of  the  jet,  is  also  unknown.   Again,  one  equation 
for  each  unknown  is  required. 

For  each  e  unknown  not  located  on  a  boundary  (nodes  8-11, 

2 

14-17,  20-23,  26-29  in  example),  the  equation   V  8  =  0   is 

applied  at  the  node  at  which  the  unknown  is  located.   The 
finite  difference  operator  shown  in  Figure  8  is  used. 

For  e  unknowns  located  at  nodes  on  the  channel  bottom 
(nodes  19-25  in  example) ,  an  equation  is  obtained  by  applying 
the  boundary  condition  —  =0   at  the  node.   The  finite 
difference  operator  shown  in  Figure  10  is  used. 

For  9  unknowns  located  at  the  channel  surface,  (nodes 
12,  18,  24,  30  in  example),  an  equation  is  obtained  by 
applying  the  boundary  conditions   t-tt  +  tt—  6  =  0   at  the  node. 

ox     K 

The  finite  difference  operator  shown  in  Figure  11  is  used. 
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The  result  is  an  equation  for  each  6  unknown.   This  system 
of  equations  is  solved  using  Gaussian  elimination  methods. 

The  method  of  solving  equations  (20)  and  (21)  is  presented 

in  the  form  of  a  flow  chart  in  Figure  15.   The  superscripts, 

p  th 

e.g.,  $*   refer  to  the  solutions  at  every  node  on  the  £ 

iteration  in  a  particular  loop.   The  quantity  Rq(^*  ,8  )  is 


defined  as  the  right  hand  side  of  equation  (20)  based  on  the 


m       n 
known  values  of  \\>*      and  6  .   The  quantity  R   is  the  right  hand 


side  of  equation  (21)  defined  similarly. 

The  convergence  tests  for  iJj*  and  e  are: 


i-1 


e1"1  -  e1 


e 


l-l 


<   0.001 


<   0.001 


It  is  assumed  that  when  the  above  criteria  are  satisfied 
at  every  node,  the  solution  has  converged. 

It  should  be  pointed  out  that  the  governing  equation  is 
not  applied  at  every  unknown  node.   The  non-linear  terms 
calculated  from  the  values  of  functions  from  the  previous 
iteration  affect  the  known  value  in  the  governing  equation  only 
at  nodes  at  which  the  governing  equation  is  applied.   However, 
the  solution  for  every  node  is  affected,  because  the  equations 
are  solved  simultaneously  by  Gaussian  elimination  techniques. 

In  the  calculation  of  the  right  hand  sides  of  equations  (20) 
and  (21) ,  two  different  techniques  were  used.   The  first 
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Figure  15  -  Flow  chart  for  iterative  scheme 
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technique  utilized  the  operators  shown  in  Figure  6,  7,  and  9, 
which  were  derived  earlier. 

A  second  method  ,  which  provided  improved  numerical 
stability,  was  also  used.   This  method  is  described  in  [13]-. 
This  is  simply  an  alternate  way  of  calculating  the  Jacobian, 
which  appears  both  in  the  right  hand  side  of  equation  (20) 
and  equation  (21) . 

Using  this  method,  the  Jacobian  Operator  is  expressed  in 
three  alternate  forms: 


J(A'a)  9X  3Y    3Y  3X 


=  i_  tA  i£)  -  i_  (A  i!) 

3X  K       3Y'    3Y  v   3X; 


3Y  y       3X;     3X  V   3Y; 


When  the  quantities  are  grouped  in  this  manner,  applying 
the  three  point  central  difference  operator  shown  in  Figure  9 
leads  to  three  different  numerical  methods  of  calculating 
J(A,B).   By  calculating  J(A,B)  by  each  of  these  methods,  and 
using  the  average  of  the  three  values,  some  forms  of  numerical 
instability  are  eliminated. 

This  method  may  be  applied  directly  to  equation  (20) . 

However,  in  order  to  apply  it  to  equation  (21) ,  the  value  of 

2 

V  ijj*  must  be  calculated  at  each  node.   This  is  done  by  applying 

the  finite  difference  operator  shown  in  Figure  8.   The  method 

2 

described  above  can  then  be  applied  to  calculate  J  ( ^ * ,  V  \p*)  . 

The  computer  program  used  to  obatin  the  results  which  are 
presented  is  included  in  Appendix  B. 
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IV.   RESULTS 

A.   PRESENTATION  OF  RESULTS 

In  the  results,  velocity  data  are  presented  in  the  form 
of  streamlines  in  the  channel.   Temperature  data  are  presented 
in  the  form  of  isotherms  in  the  channel.   Both  the  streamlines 
and  the  isotherms  are  presented  in  non-dimensional  form,  and 
are  plotted  with  respect  to  non-dimensional  values  X  and  Y. 
These  graphs  are  based  on  linear  interpolation  between  calcu- 
lated values  of  ty*    and  6  at  nodes  in  the  grid.   It  should  be 
noted  that  in  the  graphical  plots  for  both  ty*   and  6,  the  scale 
in  the  Y-direction  is  twice  that  in  the  X-direction.   The  grid 
used  for  all  calculations  consisted  of  7  nodes  in  the  Y- 
direction  and  19  in  the  X-direction. 

For  all  results  presented,  the  surface  boundary  condition 
on  the  temperature  is  —  +  t—   6  =   0  .   Additionally,  the 
ratio  of  the  flow  rate  of  the  jet  to  the  flow  rate  of  the 
river  is  0.2  in  all  cases.   This  value  was  chosen  as  a  repre- 
sentative value  for  an  actual  power  plant,  and  was  felt  to  be 
small  enough  so  that  the  surface  boundary  condition  on  the 
vertical  velocity  (i.e.,  no  waves)  remains  valid. 

In  Figures  16-21,  the  streamlines  and  isotherms  for  a 
transverse  jet  of  width  5D/6  are  presented.   The  Reynolds 
number  and  Prandtl  for  all  these  figures  are  Re  =  0.3  and 
Pr  =  6.0.   Results  for  three  different  values  of  the  ratio 
of  the  Grashof  number  to  the  Reynolds  number  are  presented; 
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Gr/Re  =  600,  Gr/Re  =  0,  and  Gr/Re  =  -300.   The  case  of  a 
negative  value  of  Gr/Re  implies  that  the  temperature  of  the 
jet  is  less  than  the  temperature  of  the  stream.   Hence,  the 
buoyant  force  acts  in  the  opposite  direction,  i.e.,  downward. 

In  Figures  22-27,  the  streamlines  and  isotherms  for  a 
longitudinal  jet  of  width  of  D/6  located  between  Y  =  1/3 
and  Y  =  1/2  are  presented.   For  all  these  figures,  Re  =  0.3 
and  Pr  =  6.0.   The  three  values  of  Gr/Re  presented  are 
Gr/Re  =  600,  Gr/Re  =  0,  and  Gr/Re  =  -300. 

In  Figures  28-29,  the  streamlines  for  both  a  transverse 
jet  of  width  5D/6  and  a  longitudinal  jet  of  width  D/6 
between  Y  =  1/3  and  Y  =  1/2  are  presented.   The  values  of 
the  non-dimensional  parameters  for  these  figures  Re  =  40, 
Pr  =  6.0,  and  Gr/Re  =0.   No  isotherms  are  presented  because 
the  jet  is  assumed  to  be  at  the  same  temperature  as  the 
river. 

In  Figure  30  and  Table  1,  the  surface  temperature  as  a 
function  of  X  is  presented.   The  temperatures  for  the  longi- 
tudinal jet  are  presented  in  tabular  form  becasue  differences 
in  the  temperatures  for  the  three  different  values  of  Gr/Re 
are  difficult  to  distinguish  graphically.   For  all  cases  in 
Figure  3  0  and  Table  1  the  Reynolds  number  and  the  Prandtl 
number  are  Re  =  0.3  and  Pr  =  6.0.   The  value  of  Gr/Re  for 
each  case  is  indicated  in  the  results. 

B.   DISCUSSION  OF  RESULTS 

The  effect  of  the  buoyant  force  on  the  flow  patterns  and 
isotherms  can  be  seen  by  examining  the  graphical  results  in 
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Gr/Re  =  600.0 


0.3- 


0.2  _ 


Surface  temperature  vs.  X  for  a  transverse  jet 
of  width  5D/6,  Re  =  0.3,  Pr  =  6.0. 
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TABLE  I 


Surface  Temperature  (9  ) 
Longitudinal  Jet 


X 

V 

6 

s 

6s 

Gr/Re  =  600 

Gr/Re  =  0 

Gr/Re  =  -300 

0.0 

0.0 

0.0 

0.0 

0.833 

0.1083 

0.1077 

0.1074 

L.667 

0.0907 

0.0901 

0.0898 

2:.  500 

0.0632 

0.0635 

0.0636 

3.333 

0.0431 

0.0437 

0.0440 

4.167 

0.0294 

0.0300 

0.0304 

5.000 

0.0200 

0.0300 

0.0304 

5:.  833 

0.0137 

0.0141 

0.0144 

6.667 

0.0094 

0.0097 

0.0099 

7.500 

0.0064 

0.0066 

0.0067 

8.333 

0.0044 

0.0046 

0.0047 

9.167 

0.0030 

0.0031 

0.0032 

10.000 

0.0021 

0.0022 

0.0022 
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Ficprares  16-21  and:  22&-2J~,..     TJies  case  of  Gr/Ree==0C implies  the 
buoyant  effect  is  negdJigdiiler  iir.  the:  momentum.- equation;   In 
tlnis  case,,  ncc  biinyantt  f  rrroBs  ac±s.=  on  thes  fibwy  ,  aadi  the  :momen turn 
equation,  is  independent  off  thee  temperatures  equation  ^ . ,  A  posi- 
tive value  of  Qi/Rfe  impdliiess  that- the  temperature eofi the  jet  is 
greater  than  T  ..   Dr.  this:  eraser  the  buoyant--  force ;aets  in  the 
positive  y~<Tt'rT3TTtri''an-<.  Sk  negative--  value: of iGf/Reeimplies  that 
the  temperature  off  the:  jett  rss  less  thanr.Tl^.  .  Inr.this  .-case, 
the  buoyant:  force  acts?  iir  thee  negative  Yrdirection;. .  In  order 
to  examine  the  efTeact:  off  buoyancy ;,  the:  casessihrwhich  Gr/Re 
is  men— zero  are  compared:"  toe  the:  case  in  which  .Gr/Re  =  0 . 

For  the  transverse:  jet,,  the  ty*   —  0  streamline : indicates 
the  upper  boundary  afr  the  flow  introduced:  by  the :  jet.   No 
flaw  ever  crosses  this  boundary..  Because  the :f low  introduced 
by  the  jet  occupies  the  space:  below  this  boundary,  the  flow 
originally  in.  the  channel:  iss  flowing  through:. a:. smaller  area, 
and!  is  therefore  accelerated:  downstream  of £  the: jet z  . 

A  positive  value  of"  Gr/Re causes  the;  transverse  jet  to 
penetrate  further  into  the  stream  than  the_ j et -for : which 
Gr/Re  =  Q,  particularly  near  the  jet  location..  As -a  result 
of  this  larger  "obstruction"  to  the  flow  in  the  channel,  local 
regions  of  decelerated  and  accelerated  flow  are  created  near 
the  jet  (Fig.  17).   A.  negative  value  of  Gr/Re  restricts  the 
penetration  of  tire  jet:..  Dr  this  case,  the:  flow  originally 
im   the  channel  is  able  tor  surpass  the_  jet- inr a: smooth  manner 
(Fig.  13)) .. 
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Fear  tliies  Bxn^i±lid±n:all  jet^  ,  thee  ef  feet -.of  the  buoyant 
force  car.  the  ffLawis:  not-  ass  pronounced  ."as  in  the  transverse 
Jet..  Th±s?  rsE  because  the-  initial! velocity  of  the  jet  is  in 
the?  dawnsstreairr  direction,  audi  therefore^  the  heat  added  by 
the  jet:  is£  corrverrted.  downstream-  more  -quickly  than  for  the 
transverse:  jeefc..  Hence,  the:  local- values  of  8  are  not  as 
great:  irr.  the:  case:  of:  the'  longitudinal  1  jet.   Since  the 
buoyant:  fame:  is:  proportionall  tor  6  ,  .  the-:buoyant  forces  near 
the  jet  SB2E£  ncct  as:  great  in-  the: case: of  a  longitudinal  jet. 

Tlie  effect  ocf  the  buoyant:f brce:can:  be  seen  by  comparing 
the  streamlines:  fox  the  case:  in:  which  Gr/Re  is  greater  than 
zero  to  tire  case:  in  which  Gr/Rer—O.  .  The  positive  buoyant 
farce  causes-  the  streamlines  to- be' bent  upward  with  respect 
ta  the  case:  in.  which  Gr/Re  —0'  (Fig;  .22).   This  is  due  to 
the  velocity  in  the  positive'  Y'-direction  caused  by  this 
buoyant:  farce..  A\  negative  value  of  :"G£/Re  has  the  opposite 
effect  (F"ig\.  231),.  causing  the:  streamlines  to  be  bent 
down-ward.. 

The  effect:  of  buoyancy  on  the  temperature  field  for 
both  jets  can  be  seen  by  examining  the  isotherms  presented 
in  Figures  19—21  and  25-27.   At  downstream  locations  near 
the  jet,  temperatures  are  greater  for  a  positive  buoyant 
force  than  for  the  case  in  which  there  is  no  buoyant  force 
present..  Th±s  is  due  to  the  fact  that  a  positive  buoyant 
farce  contributes:  to  the  heat:  transfer : in  the  vertical 
<ff  rectiair  by  increasing  the:  vertical  Lvelocity  component. 
Hence,  the  heat  added  by  the  jet  is  dispersed  more  quickly 
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£od  tire  vextrcall  direction, .  andi  temperatures  near  the  jet 
ai^  frxgherr..  Tlriiss  results;  inr.anr.  increased  surface  temperature 
Because  aTX.  heat:  added  by the;  jet-  must  leave  the  channel 
ttirrmigii.  the  surfaace^  and;  the;  amount:  of  heat  transferred  at 
the:  surface:  rs  propoxtionall  to:  the: surf ace  temperature,  more 
heat  is  txarrsferrecL  out  off  the: channel  near  the  jet  for  a 
pxssitiLve  buoyant  force  thanr.fDrra;  zero  buoyant  force. 
TJJrerjgfore  thee-  temperatures;  ait  downstream  locations  are  lower 
for:  the  positive:  buoyant  force,.  This  effect  can  be  seen  by 
examining  the:  surface  temperature : data  presented  in  Figure 
3E0..  This  effect:  occurs  for  both,  the  longitudinal  and  trans- 
verse: jets,,  but:  is~  more  pronounced  for  the  transverse  jet. 
R.   negative  buoyant  force  has  the  opposite  effect,  causing 
tire  temperatures"  near  the  jet  to:  be  smaller ,  and  the  down- 
stream tempexatnres  higher. 
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w..    discussion' 

A.  CQmEEBCBSSGK   0E"  ETNEKR-  SOLUTION. '  FDRr  STREAM  FUNCTION 

The  validity  of  thee  raatxix  of f  llnearr equations  for  the 
stream  function  can  be verified  by;  constructing  a  conver- 
gence plat  to:  determines  the  value  to"  which  the  numerical 
solution  would:  crmveacge:  iff  an  infinite:  number  of  nodes  were 
raised.   This7  value  may  then  be  compared:  to:  an  exact  solution. 
2l  canv engence  plot:  is?  constructed: byy plotting  the  value 

of  the  numerical  solution  obtained  at*:  a:  given  point  in  the 

2 

field  of  the  problem  versus  1/N  ,  where;  N  is  the  number  of 

divisions  per  side  af:  the  field.   By  making  the  number  of 
divisions  greater,,  different  values:  of  -  "the  numerical  solution 
are  obtained  at  the  same  point  in  the -field.   The  plot  of 
the  value  of  the  solution  versus  1/N2- may  then  be  extrapo- 
lated to  determine  the  point  at  which,  the.  plot  .crosses  the 
1/N2  =  0  axis..   This  value  of  the  solution  is  then  taken  as 
the  value  to  which,  the  numerical  solution  would  converge  if 
an  infinite  number  of  divisions  were  possible. 

In  order  to  apply  this  method  to  the  linear  ( zeroth 
order)  solution  for  the  stream  function,  a  fixed  field  must 
be  chosen.   The  field  chosen  was  a  section  of  the  channel 
of  length  5D.   The  case  considered  was  a  longitudinal  jet 
of  width  2H  located  at  the  upstream  boundary  and  in  the 
center  of  the  channel...  By-  utilizing  a:  grid  in  which  the 
spacing  in  the  X'-direction  is  five-  times  the  spacing  in  the 
Y-direction,  a  grid  with  an  equal  number  of  divisions  on 
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each  edge  of  thes  Sjedxt  isi  possibles  .  CSnvergence  plots  were 
constructed  at  the:  points:  X  '=  2*  5^  ,  YY==(K  25  and  X  =  2.5, 
Y  =  Q.T5.  The  ntesh:  sizes'  that-  wereeused: were  N  =  4,  N  =  8, 
and  K  =  12..  The?  N"  =  44  grid  was-  too:  coarse  to  yield  useful 
convergence  data..   Since  a  grid:  corresponding  to  N  =  16 
requires  excessive  computer-  storage.,  .  only  the  solutions 
<narxespo:uding;  to:  N:  =  8T  and"  N  =-- 12' were,  used  in  the  extrapo- 
lation process..  The  convergence- plots: for  these  two  points 
axe  shown  in.  Eignres:  32.   and  32.'. , 

An  exact  sccluticcn  at  these:  two:  points  is  obtained  by 
noting  that  the  surface  and  bottonr  boundary  conditions  for 
the  velocity  are  the  same  as  those:  for  Couette  flow  between 
flat  plates..   The  specification  of :  a:  known  flew  rate  for 
incompressible  flaw  is  equivalent  to: specifying  a  driving 
pressure  gradient  in  the  downstream  direction.   Hence,  for 
flow  in  the  channel!  far  enough  downstream  so  that  the  effect 
of  the  jet  on  the  flow  has  been  damped  out,  the  exact 
solution  is  the  solution  for  Couette  flow  between  flat 
plates  with  an  imposed  pressure  gradient.   This  solution  is 
given  in  Schlicting  [14],   It  should  be  noted  that  this 
velocity  profile  is  parabolic  in  y,  and  that  U  (0)  =0  and 
U  (D)  =  Us.   For  an  incompressible  flow,  the  pressure 

gradient  UL  may  be  determined  from  a.  specified  flow  rate. 
3X 
At  the  downstream  boundary,  several  assumptions  are  made. 

It  is  assumed  that  the  effect  of:"  the:  jet  may  be  neglected,  . 

with  the  exception  of  the  increased:  flow  added  by  the  jet, 

and  that  v  =  0..   It  is  also  assumed  that  the  velocity  profile 
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Figure  31  —  Convergence  plot  for  stream  function  at 
X  =  2.5,  Y  =  0.75. 
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Figure  32  —  Convergence  plot  for  stream  function 
at  X  =  0.25,  Y  =  2.5. 
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is  prarrTj-frryT; iks  rrr.  yy,,  arret  thafc.  th.ee  fibw>rateeis:known.   In 
acM£trfnnT./;  the  wrl rrrri ±y •  a±^  y  -  =•  0  is?  equal  "-to:  zero  and  that 
the  velocity:  aEt  jy  =  CD  iss  equal!  to:  U"s,  .  Since: the  velocity 
gaiCEfff  I'g-  is:  pa^a&?o3  xcr.  irr.  y>..  and'  sincee  the --three  conditions 
wliich.  deter  mine;  the:  crreffnjcients;  of :"  the-:parabola  are 
identical-  to  the:  conditions-  for  Couette  flow  with  an  imposed 
pressure:  grsdieartt, ,  the-  velocity  prof  ile-atithe  downstream 
fecrurrdary-  is:  id'enr-icarlL  tec  the  exact:  solution; 

Tfc  was:  assumed"  that:  the'  exact- solution:  was  valid  at 
X  =  2..5..  The  varlu.es:  obtained  from- the^  convergence  plot  may 
then,  be?  compared"  tor  those:  obtained  from:  the  exact  solution. 
At  the  point  X:  =  2^5r,  Y.'  =  0.25,  the  exact: solution  is 
ty*   =■   Q. .12723:;  the  solution  obtained,  from-  the  convergence 
plot  is:  ^*  =  Of.  .127331..  At:  the  point  -  X"  ==  2  .  5  ,~  Y  =  0.75,  the 
exact  solution  is-  ip*   =  0'.. 83170;  the.  solution  obtained  from 
ttfte  convergence:  ploct:  is:  ty**  -   0.83163^  .  As:a:-result  of  the 
agreement  between. these  values,  it_may  be .concluded  that 
the  linear-  solution  fear  ifr**,  which  is- the  zeroth  order  solu- 
tion for  the  iterative  scheme,  is  valid. 

B.   NUMERICAL  STABILITY  OF  ITERATIVE  SCHEME 

By  formulating  the  problem  in  terms  of  the  linear 
solutions  as  shown  in  equations  (20)  and  (21) ,  the  nonlinear 
terms  in  the  energy  equation  are  multiplied  by  the  Peclet 
number  Be,,  where  Ee  =  Re :-.  Px..   The:  nonlinear  terms  in  the 
momentum  equation:  are  multiplied  by  Re. .  However,  using  the 
iterative  procedure  described  earlier :,  the  :  numerical  solution 
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JEaiar  the  enerjjy/  equation  becomess  unstable  at  approximately 

Fe  =  2..  Siirmes  thet  energy  andi momentum-  equations  are  coupled, 

tt&e  entire  scd2iii±an  becomess  unstablesat :  Pe  =  2. 

Severall  method's:  were:  employed:"  in:. an  effort  to  eliminate 
tEtiLs  ins^taiidilxlty . .  The  method:  initially  used  to  calculate 
the  nonlinear-  terms-  was  to:  use: thee  operators  shown  in 
Figures  5,.  (i„  and:  7..  Usingr  this- method,  the  energy  equation 
became  unstable  at.  approximately ;-Pec=-  1.   When  uncoupled 
from  the  energy  equation,  the:  momentum  equation  became 
unstable  act  He  =  15. 

The  method"  of:  calculating-  the:  nonlinear  terms  in  terms 
of  the  iTacob±an.  aES-  described:  earlier:  was  then  attempted. 
Using  this  method ,  the  Peclet:  number:  at  which  instability 
©Grcurred  in.  the  energy  equation  was:  raised  from  1  to  2 . 
The  Reynolds  number  at  which  instability  occurred  in  the 
momentum  equation,,  uncoupled:,  from  the:  energy  equation,  was 
raised",  from  15  to  5  0. 

Decreasing-  the  grid  spacing  was : also  attempted.   This 
had  no  noticeable  effect  on  the  stability  of  the  system  of 
equations. 

An  attempt  was  made  to  begin  with  the  stable  solution 
for  Pe  =  2.0,  and  increase  the  Peclet  number  in  steps  of 
1/25.   Using  this  procedure,  instability  occurred  at  Pe  =  2.15 

C   CONCLUSIONS 

Because-  all,  solutions  presented  in  this  study  are  based 
cm  iterations-  from  the  linear:, ,  uncoupled  solutions  of  the 
derived  governing  equations,  the  validity  of  the  numerical 
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mffittiiEodEsE  usacE  toe  obtain  thesee solutions  is  of  utmost  impor- 
Hamce?..  The?.  valTidity  of::  thee  operators  rand  method  of  solution 
fixsar  trite?  Limsarr  errer.gy  eqnationr. has  5 been  demonstrated  by  its 
fwquenti  usee  with  success:-  inr.  tbei literature  [12,  15].   The 
re?sarl~ts=  off  the?  convergence'  plot-  indicate  that  the  linear 
solution  off  the?  linear  momentum'  equation  is  also  valid, 
^esfeid  Q2T.  theses  observations^  ,  it":  is  £  felt  that  the  solutions 
gareseiitlejd  rxr.  t±ds;  report:  aree  calculated  from  an  essentially 
SdOUird  basis?.. 

The?  effect  of"  buoyancy  on:,  the:  flow  patterns  and  tempera- 
ture d±str±buti on s  has  beenr  demonstrated  for  the  cases 
ccorsixfered'..  However,  the:  original -motivation  for  this 
effort,  was  tec  demonstrate  the  effect  _of  buoyancy  on  the  flow 
in?  an.  actual  physical  channels   For:  such  channels,  the 

actual  values:  of:  Reynolds  and  Peclet  numbers  experienced 

3 
a^re?  greater:  tharr  10  .   Because,  of :  the  "restraints  on  the 

range  of:  Re:  and:.  Pe.  due  to  nonlinear:  numerical  instability, 
it  was  net  possible  to  generate:  solutions  which  are  appli- 
cable to  actual  channels. 

The  maximum  value  of  the  Peclet  number  for  which  results 
were  obtained  was  approximately  2.0.   Using  the  value 
Pr  =  6  as  a  typical  value  of  the  Prandtl  number  for  water, 
the  maximum  Reynolds  number  for  which  results  for  the 
coupled  temperature  and  flow  fields  were  obtained  was 
appxcximateXy  OCX.. 

5s  a.  result  of  these  limitations;  some  of  the  assump- 
tions made  for  the  boundary  conditions  for  the  longitudinal 
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jet  may  be  questioned:..  Irr.  the  tbrmulationr.of  f  the  boundary 
condition  for:  the  temperatu  re  for-.- thee  longitudinal  jet,  it 
was  assumed:  that  nn  heat  was:  conducted:  upstream- to  the 
X-location  of:  the  jet..  This  assumption",  iss valid  for 
Pe  >>  1«   Since  the-  maximum  Peclet  nnmberr used: in  the  cou- 
pled problem  was?  Ee  =  31..8I,,  this  assumption:  is£questionable. 

For  the  velocity  bamniary  condition,  forr the : longitudinal 
jet,  it  was  assumed"  that:  the  distruhancee toe  the: flow  field 
caused  by  the  jet  does  not:  propagate- upstream'  to  the 
X-location  of  the  jet..  This;  assumption:,  issvalid  for 
Re  >>  L.   Because  the  maximum  Reynolds:. number:  used  in  the 
solution  of  the  coupled-  energy  and  momentum' equations  was 
0.3,  this  assumption  is  ncc  longer  valid. 

D.   RECOMMENDATIONS 

Using  a  finite  difference  method? of f  solution,  one 
algebraic  equation,  is  generated  for  each  node . .  In  order  to 
solve  such  a  system  of  linear  algebraic  equations,  two 
basic  methods  exist..  The  first  method  is:  the: direct  method, 
which  is  to  solve  the  simultaneous  equations  directly, 
generating  an  exact  solution  to  the  system  of  equations. 
The  second  method  is  the  indirect  method,  which  arrives  at 
a  solution  by  successive  improvements  of  an  approximate 
solution.   This  process  is  repeated  until  a  desired  conver- 
gence criteria  is  satisfied. 

The  basicr  method  off  soiution  used."  in  this  study  was  the 
direct  method..   Using-  this:  method,'  the  number:  of  calcula- 
tions required  to  achieve  a  solution  is  known,  which  is  not 
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true  erf  the  indirect:  method:. .  The-  directr method : is  also  well 
adapted  to  machine  aaelcxuLacticms  .   Ihr.  additibnr,  the  Gaussian 
elimination:  techniques  of:  scdJution  by; "thee direct: method  is 
the  most  efficient:  method:  erf'  solvingr  a?.  systemrof  simultane- 
ous Linear  aigehr.aic:  equations.   Usingr thisi method,  however, 
extensive  campu-terr  storage  is  required..  For:-  example ,  a 
finite  difference  solution,  for  a  gride  involving  10  nodes  in 
each  direction-  generates;  10D  simultaneous:  equations  for 
each  unknown...   In.  order-  trr  solve  thisz  prDblerarby  the  direct 
method,,  a  1GCF  by:  IOlOC  matrix  must  be:  stored  for:  each  unknown. 

On.  the  basis  of  this  work ,  the:  solution:  of :  nonlinear 
algebraic  equations  using-  the  iterative:  method  proposed 
earlier  is  not  feasible  when  the  value- of :"  the  nonlinear 
terms  becomes  large..  For- the  solution: of :  such  a  problem, 
the  indirect  method  of  solution  is  recommended. 

One  passible  method"  af'  solving^  such-  a:  system  of  nonlinear 
algebraic  equations  is-  to  employ  a  numerical'  relaxation 
scheme.   This  method  is  discussed  in  [12]_.   One  advantage 
of  this  method  is  that  once  the  relaxation  pattern  is 
established,  the  need  for  storing  matrices  is  eliminated. 

Another  possible  method  of  solution  is  to  recast  the 
problem  in  the  form  of  vorticity  transport 

0)   =   V2i|;*  (22) 

Ere  t*£l  ^  -  *£l!f!j  =   V2W.  -££lAL  (23) 

2T  ax:   ax:  8Y  Re;  3Y: 

Re  Pr   [1*1  IL-Ml.il]   =   V2e  (24) 

3Y   3X    3X   3Y 
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where  cu  is  the  varticitv..  Thee  finite:  ddfferencec  forms ;  of : 
these  equations  can.  fee  scllved:  fim' tiiee  vadhee- off  thee  unknown 
at  a  node  in  terms  of  the  vaXueas  o±f  o±h©rr  unknowns  rat  -_ 
surrounding  nodes-   Fear  instance:,,  equation:  (22J,  can r. bee 
solved  for  cu  in  terms  off  ^  at:  surrxrunding-  points;  .  The: 
values  obtained  from  the  previous  iteration-,  canr.bef  used 
to  calculate  new  values..   B:y •  chonsingr  some--  appropriate  _ 
starting  point  p  an  iterative:  solutirxn.  scheme:,  canr. bee 
established. 

A  further  recommendation,  is  that:  aa  gradedc  network  be 
used  in  the  numerical  solntion..  This:  graded-  network  should 
have  finer  divisions  near'  the  jet-  location. .  Thisrwould: 
improve  the  accuracy  of  the  method  of'  modeling  the  jet, 
and  also  provide  more  grid"  points  near:  the-  jet>  .  where  the 
maximum  changes  in  fTaw  and  temperature:  are:  experienced. 
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APPENDIX  A 
COMPUTER  PROGRAM  USED  TO  DEVELOP 
FINITE  DIFFERENCE  OPERATORS 
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APPENDIX' :  B3 

GHMEunnr  program  used~  in:  numerical  solution 
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